rm(list=ls())

library(tidyr)
library(tibble)
library(stargazer)
library(jtools)
library(dplyr)
library(haven)
library(readstata13)
library(ggplot2)
library(interactions)
library(ggpubr)
library(ggplotify)

#########################################################################
# REPLICATION FILES: MILITARIZATION AND PERCEPTIONS OF LAW ENFORCEMENT  #
# JOURNAL: BJPS                                                         #
# FIGURES IN SUPP. APP. FOR:                                            #
# ADDITIONAL CONTINUOUS INTERACTION EFFECTS                             #
#########################################################################

#### STEP 1: LOAD RESULTS                           ####


#### STEP 1: READ DATA                              ####

militconj <- read.dta13("MilitConjoint.dta")
baselines <- list()

#### STEP 2: DECLARE THEME                          ####

theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 11, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=0.5), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=1.5),
      axis.title.x =      element_text(size =9),
      panel.background = element_rect(fill = "#f7f7f7"),
      strip.background = element_rect(colour="black", fill="white"),
      legend.position="bottom"
    )
}

#### STEP 3: FIGURE A.23: INTERACTION BETWEEN MILITARY PRESENCE AND ASSAULT WEAPON ATTRIBUTE   ####

# EFFECTIVENESS
lm_effev <- lm(zeffective ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jev_weap<- johnson_neyman(lm_effev, 
                          pred = weapon, 
                          modx = countevents, 
                          alpha = .05,
                          insig.color="grey",
                          sig.color="black",
                          title="Effectivenesss",
                          mod.range=c(0,101))



jev_weap <- jev_weap$plot + xlab("Confrontations") + 
  ylab("Assault rifle slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)
plot(jev_weap)

jev_weap <- ggplot_build(jev_weap)
jev_weap$data[[8]]$alpha <- 0
jev_weap <- ggplot_gtable(jev_weap)
jev_weap  <- as.ggplot(jev_weap)
plot(jev_weap)

# CIVIL LIBERTIES
lm_libev <- lm(zlib ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jlib_weap<- johnson_neyman(lm_libev, 
                           pred = weapon, 
                           modx = countevents, 
                           alpha = .05,
                           insig.color="grey",
                           sig.color="black",
                           title="Civil Liberties",
                           mod.range=c(0,101))

jlib_weap <- jlib_weap$plot + 
  xlab("Confrontations") + 
  ylab("Assault rifle slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jlib_weap <- ggplot_build(jlib_weap)
jlib_weap$data[[8]]$alpha <- 0
jlib_weap <- ggplot_gtable(jlib_weap)
jlib_weap  <- as.ggplot(jlib_weap)

plot(jlib_weap)


# CORRUPTION
lm_correv <- lm(zcorrupt ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)

jcorr_weap<- johnson_neyman(lm_correv, 
                            pred = weapon, 
                            modx = countevents, 
                            alpha = .05,
                            insig.color="grey",
                            sig.color="black",
                            title="Corruption",
                            mod.range=c(0,101))

jcorr_weap <- jcorr_weap$plot + 
  xlab("Confrontations") + 
  ylab("Assault rifle slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jcorr_weap)


# NEIGHBORHOOD
lm_neighev <- lm(zneighb ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jneigh_weap<- johnson_neyman(lm_neighev, 
                             pred = weapon, 
                             modx = countevents, 
                             alpha = .05,
                             insig.color="grey",
                             sig.color="black",
                             title="Neighborhood",
                             mod.range=c(0,101))

jneigh_weap <- jneigh_weap$plot + 
  xlab("Confrontations") + 
  ylab("Assault rifle slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jneigh_weap <- ggplot_build(jneigh_weap)
jneigh_weap$data[[8]]$alpha <- 0
jneigh_weap$data[[9]]$alpha <- 0
jneigh_weap <- ggplot_gtable(jneigh_weap)
jneigh_weap  <- as.ggplot(jneigh_weap)

plot(jneigh_weap)

# JOIN

pres_jn_weap <-  ggarrange(
  jev_weap, jlib_weap, jcorr_weap, jneigh_weap,
  common.legend = TRUE, legend = "bottom", ncol = 4
)

plot(pres_jn_weap)

#ggexport(pres_jn_weap, 
         #filename="Presence_JNeum_Weapon.png",
         #width = 1400,
         #height = 400)


#### STEP 4: FIGURE A.24: INTERACTION BETWEEN MILITARY PRESENCE AND MALE ATTRIBUTE   ####

# EFFECTIVENESS
lm_effev <- lm(zeffective ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jev_gen<- johnson_neyman(lm_effev, 
                         pred = gender, 
                         modx = countevents, 
                         alpha = .05,
                         insig.color="grey",
                         sig.color="black",
                         title="Effectivenesss",
                         mod.range=c(0,101))

jev_gen <- jev_gen$plot + xlab("Confrontations") + 
  ylab("Male slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)
plot(jev_gen)

# CIVIL LIBERTIES
lm_libev <- lm(zlib ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jlib_gen<- johnson_neyman(lm_libev, 
                          pred = gender, 
                          modx = countevents, 
                          alpha = .05,
                          insig.color="grey",
                          sig.color="black",
                          title="Civil Liberties",
                          mod.range=c(0,101))

jlib_gen <- jlib_gen$plot + 
  xlab("Confrontations") + 
  ylab("Male slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jlib_gen <- ggplot_build(jlib_gen)
jlib_gen$data[[8]]$alpha <- 0
jlib_gen$data[[9]]$alpha <- 0
jlib_gen <- ggplot_gtable(jlib_gen)
jlib_gen  <- as.ggplot(jlib_gen)

plot(jlib_gen)


# CORRUPTION
lm_correv <- lm(zcorrupt ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)

jcorr_gen<- johnson_neyman(lm_correv, 
                           pred = gender, 
                           modx = countevents, 
                           alpha = .05,
                           insig.color="grey",
                           sig.color="black",
                           title="Corruption",
                           mod.range=c(0,101))

jcorr_gen <- jcorr_gen$plot + 
  xlab("Confrontations") + 
  ylab("Male slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jcorr_gen <- ggplot_build(jcorr_gen)
jcorr_gen$data[[8]]$alpha <- 0
jcorr_gen <- ggplot_gtable(jcorr_gen)
jcorr_gen  <- as.ggplot(jcorr_gen)

plot(jcorr_gen)


# NEIGHBORHOOD
lm_neighev <- lm(zneighb ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jneigh_gen<- johnson_neyman(lm_neighev, 
                            pred = gender, 
                            modx = countevents, 
                            alpha = .05,
                            insig.color="grey",
                            sig.color="black",
                            title="Neighborhood",
                            mod.range=c(0,101))

jneigh_gen <- jneigh_gen$plot + 
  xlab("Confrontations") + 
  ylab("Male slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jneigh_gen)

# JOIN

pres_jn_gen <-  ggarrange(
  jev_gen, jlib_gen, jcorr_gen, jneigh_gen,
  common.legend = TRUE, legend = "bottom", ncol = 4
)

plot(pres_jn_gen)

#ggexport(pres_jn_gen, 
         #filename="Presence_JNeum_Gender.png",
         #width = 1400,
         #height = 400)

#### STEP 5: FIGURE A.25: INTERACTION BETWEEN MILITARY PRESENCE AND LIGH SKIN ATTRIBUTE   ####

# EFFECTIVENESS
lm_effev <- lm(zeffective ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jev_skin<- johnson_neyman(lm_effev, 
                          pred = race, 
                          modx = countevents, 
                          alpha = .05,
                          insig.color="grey",
                          sig.color="black",
                          title="Effectivenesss",
                          mod.range=c(0,101))

jev_skin <- jev_skin$plot + xlab("Confrontations") + 
  ylab("Light skin slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)
plot(jev_skin)

# CIVIL LIBERTIES
lm_libev <- lm(zlib ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jlib_skin<- johnson_neyman(lm_libev, 
                           pred = race, 
                           modx = countevents, 
                           alpha = .05,
                           insig.color="grey",
                           sig.color="black",
                           title="Civil Liberties",
                           mod.range=c(0,101))

jlib_skin <- jlib_skin$plot + 
  xlab("Confrontations") + 
  ylab("Light skin slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jlib_skin)


# CORRUPTION
lm_correv <- lm(zcorrupt ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)

jcorr_skin<- johnson_neyman(lm_correv, 
                            pred = race, 
                            modx = countevents, 
                            alpha = .05,
                            insig.color="grey",
                            sig.color="black",
                            title="Corruption",
                            mod.range=c(0,101))

jcorr_skin <- jcorr_skin$plot + 
  xlab("Confrontations") + 
  ylab("Light skin slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jcorr_skin)


# NEIGHBORHOOD
lm_neighev <- lm(zneighb ~ countevents*uniform + countevents*weapon + countevents*gender + countevents*race, data = militconj)
jneigh_skin<- johnson_neyman(lm_neighev, 
                             pred = race, 
                             modx = countevents, 
                             alpha = .05,
                             insig.color="grey",
                             sig.color="black",
                             title="Neighborhood",
                             mod.range=c(0,101))

jneigh_skin <- jneigh_skin$plot + 
  xlab("Confrontations") + 
  ylab("Light skin slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jneigh_skin)


# JOIN

pres_jn_skin <-  ggarrange(
  jev_skin, jlib_skin, jcorr_skin, jneigh_skin,
  common.legend = TRUE, legend = "bottom", ncol = 4
)

plot(pres_jn_skin)

#ggexport(pres_jn_skin, 
         #filename="Presence_JNeum_Skin.png",
         #width = 1400,
         #height = 400)

#### STEP 6: FIGURE A.26: INTERACTION BETWEEN HOMICIDE RATE AND ASSAULT WEAPON ATTRIBUTE   ####

# EFFECTIVENESS
lm_effhom <- lm(zeffective ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)
jeffhom_weap<- johnson_neyman(lm_effhom, 
                              pred = weapon, 
                              modx = thom, 
                              alpha = .05,
                              insig.color="grey",
                              sig.color="black",
                              title="Effectivenesss",
                              mod.range=c(0,133))

jeffhom_weap <- jeffhom_weap$plot + xlab("Average Homicide Rate") + 
  ylab("Assault rifle slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jeffhom_weap <- ggplot_build(jeffhom_weap)
jeffhom_weap$data[[8]]$alpha <- 0
jeffhom_weap <- ggplot_gtable(jeffhom_weap)
jeffhom_weap  <- as.ggplot(jeffhom_weap)

plot(jeffhom_weap)


# CIVIL LIBERTIES
lm_libhom <- lm(zlib ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)
jlibhom_weap<- johnson_neyman(lm_libhom, 
                              pred = weapon, 
                              modx = thom, 
                              alpha = .05,
                              insig.color="grey",
                              sig.color="black",
                              title="Civil Liberties",
                              mod.range=c(0,133))

jlibhom_weap <- jlibhom_weap$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Assault rifle slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jlibhom_weap <- ggplot_build(jlibhom_weap)
jlibhom_weap$data[[8]]$alpha <- 0
jlibhom_weap <- ggplot_gtable(jlibhom_weap)
jlibhom_weap  <- as.ggplot(jlibhom_weap)

plot(jlibhom_weap)


# CORRUPTION
lm_corrhom <- lm(zcorrupt ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)

jcorrhom_weap<- johnson_neyman(lm_corrhom, 
                               pred = weapon, 
                               modx = thom, 
                               alpha = .05,
                               insig.color="grey",
                               sig.color="black",
                               title="Corruption",
                               mod.range=c(0,133))

jcorrhom_weap <- jcorrhom_weap$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Assault rifle slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jcorrhom_weap)

# NEIGHBORHOOD
lm_neighhom <- lm(zneighb ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)

jneighhom_weap<- johnson_neyman(lm_neighhom, 
                                pred = weapon, 
                                modx = thom, 
                                alpha = .05,
                                insig.color="grey",
                                sig.color="black",
                                title="Neighborhood",
                                mod.range=c(0,133))

jneighhom_weap <- jneighhom_weap$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Assault rifle slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jneighhom_weap <- ggplot_build(jneighhom_weap)
jneighhom_weap$data[[8]]$alpha <- 0
jneighhom_weap <- ggplot_gtable(jneighhom_weap)
jneighhom_weap  <- as.ggplot(jneighhom_weap)

plot(jneighhom_weap)


# JOIN

hom_jn_weap <-  ggarrange(
  jeffhom_weap, jlibhom_weap, jcorrhom_weap, jneighhom_weap,
  common.legend = TRUE, legend = "bottom", ncol = 4
)

plot(hom_jn_weap)

#ggexport(hom_jn_weap, 
         #filename="Homicides_JNeum_Weapon.png",
         #width = 1400,
         #height = 400)

#### STEP 7: FIGURE A.27: INTERACTION BETWEEN HOMICIDE RATE AND MALE ATTRIBUTE   ####

# EFFECTIVENESS
lm_effhom <- lm(zeffective ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)
jeffhom_gen<- johnson_neyman(lm_effhom, 
                             pred = gender, 
                             modx = thom, 
                             alpha = .05,
                             insig.color="grey",
                             sig.color="black",
                             title="Effectivenesss",
                             mod.range=c(0,133))

jeffhom_gen <- jeffhom_gen$plot + xlab("Average Homicide Rate") + 
  ylab("Male slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jeffhom_gen)


# CIVIL LIBERTIES
lm_libhom <- lm(zlib ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)
jlibhom_gen<- johnson_neyman(lm_libhom, 
                             pred = gender, 
                             modx = thom, 
                             alpha = .05,
                             insig.color="grey",
                             sig.color="black",
                             title="Civil Liberties",
                             mod.range=c(0,133))

jlibhom_gen <- jlibhom_gen$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Male slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jlibhom_gen <- ggplot_build(jlibhom_gen)
jlibhom_gen$data[[8]]$alpha <- 0
jlibhom_gen$data[[9]]$alpha <- 0
jlibhom_gen <- ggplot_gtable(jlibhom_gen)
jlibhom_gen  <- as.ggplot(jlibhom_gen)

plot(jlibhom_gen)

# CORRUPTION
lm_corrhom <- lm(zcorrupt ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)

jcorrhom_gen<- johnson_neyman(lm_corrhom, 
                              pred = gender, 
                              modx = thom, 
                              alpha = .05,
                              insig.color="grey",
                              sig.color="black",
                              title="Corruption",
                              mod.range=c(0,133))

jcorrhom_gen <- jcorrhom_gen$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Male slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jcorrhom_gen <- ggplot_build(jcorrhom_gen)
jcorrhom_gen$data[[8]]$alpha <- 0
jcorrhom_gen$data[[9]]$alpha <- 0
jcorrhom_gen <- ggplot_gtable(jcorrhom_gen)
jcorrhom_gen  <- as.ggplot(jcorrhom_gen)

plot(jcorrhom_gen)

# NEIGHBORHOOD
lm_neighhom <- lm(zneighb ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)

jneighhom_gen<- johnson_neyman(lm_neighhom, 
                               pred = gender, 
                               modx = thom, 
                               alpha = .05,
                               insig.color="grey",
                               sig.color="black",
                               title="Neighborhood",
                               mod.range=c(0,133))

jneighhom_gen <- jneighhom_gen$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Male slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jneighhom_gen)


# JOIN

hom_jn_gen <-  ggarrange(
  jeffhom_gen, jlibhom_gen, jcorrhom_gen, jneighhom_gen,
  common.legend = TRUE, legend = "bottom", ncol = 4
)

plot(hom_jn_gen)

#ggexport(hom_jn_gen, 
         #filename="Homicides_JNeum_Gender.png",
         #width = 1400,
         #height = 400)

#### STEP 8: FIGURE A.28: INTERACTION BETWEEN HOMICIDE RATE AND LIGH SKIN ATTRIBUTE   ####

# EFFECTIVENESS
lm_effhom <- lm(zeffective ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)
jeffhom_skin<- johnson_neyman(lm_effhom, 
                              pred = race, 
                              modx = thom, 
                              alpha = .05,
                              insig.color="grey",
                              sig.color="black",
                              title="Effectivenesss",
                              mod.range=c(0,133))

jeffhom_skin <- jeffhom_skin$plot + xlab("Average Homicide Rate") + 
  ylab("Light skin slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jeffhom_skin)


# CIVIL LIBERTIES
lm_libhom <- lm(zlib ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)
jlibhom_skin<- johnson_neyman(lm_libhom, 
                              pred = race, 
                              modx = thom, 
                              alpha = .05,
                              insig.color="grey",
                              sig.color="black",
                              title="Civil Liberties",
                              mod.range=c(0,133))

jlibhom_skin <- jlibhom_skin$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Light skin slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank()) +
  ylim(-1.5, 1.5)

plot(jlibhom_skin)


# CORRUPTION
lm_corrhom <- lm(zcorrupt ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)

jcorrhom_skin<- johnson_neyman(lm_corrhom, 
                               pred = race, 
                               modx = thom, 
                               alpha = .05,
                               insig.color="grey",
                               sig.color="black",
                               title="Corruption",
                               mod.range=c(0,133))

jcorrhom_skin <- jcorrhom_skin$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Light skin slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

plot(jcorrhom_skin)

# NEIGHBORHOOD
lm_neighhom <- lm(zneighb ~ thom*uniform + thom*weapon + thom*gender + thom*race, data = militconj)

jneighhom_skin<- johnson_neyman(lm_neighhom, 
                                pred = race, 
                                modx = thom, 
                                alpha = .05,
                                insig.color="grey",
                                sig.color="black",
                                title="Neighborhood",
                                mod.range=c(0,133))

jneighhom_skin <- jneighhom_skin$plot + 
  xlab("Average Homicide Rate") + 
  ylab("Light skin slope") + 
  theme_bw1() + 
  theme(legend.title = element_blank(), legend.position = "none") +
  ylim(-1.5, 1.5)

jneighhom_skin <- ggplot_build(jneighhom_skin)
jneighhom_skin$data[[8]]$alpha <- 0
jneighhom_skin <- ggplot_gtable(jneighhom_skin)
jneighhom_skin  <- as.ggplot(jneighhom_skin)

plot(jneighhom_skin)


# Join

hom_jn_skin <-  ggarrange(
  jeffhom_skin, jlibhom_skin, jcorrhom_skin, jneighhom_skin,
  common.legend = TRUE, legend = "bottom", ncol = 4
)

plot(hom_jn_skin)

#ggexport(hom_jn_skin, 
         #filename="Homicides_JNeum_Skin.png",
         #width = 1400,
         #height = 400)


